Calculating response functions in time domain with non-orthonormal basis sets 
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I. INTRODUCTION 

As first-principles calculations become more and more 
important in various research fields such as physics, 
chemistry, materials science, and recently geology and 
biology, the demand for calculation of larger and larger 
systems is growing rapidly. One of the answers to this de- 
mand is the so-called order- methods, which compute 
the electronic band structure, the total energy, and other 
quantities with computational time and storage propor- 
tional to N, the number of the atoms in the system. For 
very large systems, these methods are much faster than 
the conventional diagonalization methods, which require 
computational efforts proportional to A^'^. 

The order- A^ methods may be classified into two steps. 
The first step is minimizing the total energy to obtain 
the ground state of the self-consistent one-particle Hamil- 
tonian. The second step is extracting dynamic proper- 
ties such as linear and nonlinear-response functions from 
this Hamiltonian. While the first step has been exten- 
sively studied and also comprehensive reviews are 
available ||9|,|l0|, the second step has been studied by only 
few papers ||ll|-[l5|], including the particle source method 
p6| , p7t and the projection method [p"8[-|2l|, which use the 
numerical solution of the time-dependent Schrodinger 
equation |22|, and projected random vectors [ p3[ . 

The purpose of this Rapid Communication is to extend 
the formalism of the projection method to nonorthonor- 
mal basis sets p^-p^ , on which many order- A^ total en- 
ergy minimization methods are built, so that the full ab 
initio calculation from the total energy minimization to 
the response function is possible. 



II. NONORTHONORMAL BASIS SET 

In this section, let us review the description of a sys- 
tem with a Hilbert space spanned by finite numbers of 
linearly independent nonorthonormal bases {jc^ct)}. We 



distinguish a vector in the Hilbert space from its com- 
ponents by using the braket notation for a vector in the 
Hilbert space and the tensor notation and the matrix 
notation [ p6| for its components. 

The overlap matrix is defined as a Hermitian matrix 
with subscripts. 



Sap = iVal^p) = Sp^ 



(1) 



Then the inverse matrix is defined as a matrix with su- 
perscripts that satisfies 
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where 5" is Kroneker's delta. Then the dual basis set 



((/?" I is defined by 



(3) 



which is used only in formal description, but not in real 
numerical calculations. These two basis sets are biorthog- 
onal and bicomplete, 
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Ei^")('^"i = ^' 
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where / is the identity operator. 

An arbitrary state |(/)) can be expressed in original or 
dual basis set, 

10) = E'^"!^") = Y,<j^'^\v^)Spa = Y,Mv''). (6) 

a a.J3 

where and (pa are the components in each basis set, 
which are related to each other by 



^J2spar- 
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The components of |0) are represented by a column vec- 
tor (p = [(j)^, 0^, • • • , (/)^] , where t indicates the trans- 
pose of a vector or matrix, and its dual (0| is represented 

by a row vector (p = [0* , 02 1 ■ ' ■ j 4'*n] ■ 

The lower-indexed components of an operator, the 
Hamiltonian H for example, are defined in the original 
basis set by 
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Hal3 = {(pa\H\ip0). 



(8) 



III. RANDOM VECTORS 



Then the mixed-indexed components are defined by 



7/3- 



(9) 



The manipulation of state vectors and operators is most 
conveniently expressed in the mixed representation. For 
example, IV'} ~ H\(p) becomes — J2f3 -^"p'l^^ ■ There- 
fore, we can introduce the matrix notation, xp = Hcj) 
where the bar over the matrix symbol indicates the raise 
of the first index H = {i?"^}. Then Eq. (|) is rewritten 



H = S^H, 



(10) 



where H is the matrix {Hap}. Now H is not Hermitian 
matrix anymore, since 



= HS-i = SHS 1 ^ H. 



(11) 
(12) 



Note that the full calculation of S~^, which costs 0{N^) 
CPU time, is not necessary to obtain a good approximant 
of H from a sparse H |2^,^ . One of the advantages of 
H over H is that power of H is easily calculated without 
exphcitly multiplying S^^ pi], 



The matrix form of the eigenvalue problem 



becomes 

He/, (Ef,) = £;^0 (Ep) 
and the dual of Eq. ( [l6| ) becomes 

<P {Ef3)B. = Ep^ (Ep). 



(13) 
(14) 

(15) 

(16) 
(17) 



The eigenvectors, Eqs. ( |16| ) and (17), define the eigen- 
states 



a 

{Ep\=Y,4>a{Ep){^% 



(18) 
(19) 



which satisfy the biorthonormality and the bicomplete- 
ness 



Y,\E^){E^\^I. 



(20) 
(21) 



Let us define random states |29| , |30|] by 

P 



(22) 
(23) 



where {ItyC/j)} and {(<y5"|} are the basis set used in the 
computation and its dual basis set, respectively. 
Their components 



(24) 



are the pseudorandom numbers that satisfy the statistical 
relation 



(25) 



where (( • )) indicates the statistical average. Note that 
the transformation of the random vector to its dual does 
not contain the overlap matrix S in Eq. (^4|), unlike the 
general rule for usual vectors in Eq. (|^). 

These random vectors may be also expressed by the 
eigenstates of H by substituting Eq. ( ^ ) into Eqs. (|22" 
and (I2F 



/37 p 

{^^Y.^l{v'\E^){E^\=Y.C{E^\, 
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where 



C/3 = Y{Ep\ip^)£,^, 
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(26) 
(27) 

(28) 
(29) 



Although we do not know the actual value of Cai Cpi 
{Ep\, or \Ea): we can derive the statistical relation of the 
random variables C,p as follows: 

= Y,{EpW-,){v''\E^) - {Ep\E^) = So^p. (30) 



This relation is very important, as we will see later. 

One of the useful features of random states is that the 
expectation value of an operator X in terms of the ran- 
dom states gives trace of the operator. 



a,P a 



which is identical to the trace calculated with an or- 
thonormal basis set In) because 
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[x) ^ tT[X] ^ {n\'Pc.)X'^i3{'P^\n) (31) 
= ^(^'^|^„)X"^ = ^X«. (32) 

Q,/3 a 



IV. PROJECTED RANDOM VECTORS 

Then the projected random vectors are defined by 

=e{Ef-H)^ =J2cmi'm (33) 

m 

= * e{Ef - H) = ^ ^ (34) 

m 

where Cm are the coefficients for the Chebyshev polyno- 
mial expansion of the step function 



d{x) 



(x < 0) 

1 {x > 0). 



(35) 



The random vectors multiplied by the Chebyshev poly- 
nomial Tm{Ji) 



tA„, = T,„(H)* , 

^ „, = * r„,(H), 

are calculated by using the recursion formulas 

m+l = 2Hl/j „ - 
^ m+l = 21/^ mH-l/i 



(36) 
(37) 



(38) 
(39) 



The coefficient vectors, Eqs.(|33|) and (eq:step.def2), de- 
fine the projected random states 

l$ij,)EE^|^„}($^,)"= ^ \Ep)Cp, (40) 



E0<Ef 



One of the useful features of projected random states is 
that the expectation value of an operator X with them 
gives the trace of the operator over the Fermi occupied 
states, 



V. TIME EVOLUTION 

The time-dependent Schrodinger equations corre- 
sponding to the eigenvalue equations (^6|) and (^7|) be- 
come 



i-cf> (t)-H</. (t), 
(t) ^ (t)H. 



(44) 
(45) 



The formal solutions of the time-dependent equations be- 
come 



<f> (i) = e-^«> {t - 0), 
4> (t) = (P {t = 0)e+'^K 



(46) 
(47) 



For numerically calculating the time evolution of the co- 
efficients, we use the leap frog method [p2|, 

(f) {t + At) ^ -2iAm4> {t) + 4> (t- At), (48) 
<jy{t + At)^+2iAt4>{t)H + (^{t-At), (49) 

where At is the time step. 



VI. LINEAR RESPONSE FUNCTION 

When an impulse of perturbation AS(t) is applied to 
the system described by the Hamiltonian H, the time 
evolution of the wave function is described by the time- 
dependent Schrodinger equation in the matrix form 



it) = {H + ASit)}^ (t), 
-z^* {t)^4 {t){n + A6it)}, 



(50) 
(51) 



where A — A is the matrix of A in the mixed repre- 
sentation. Note that the impulse AS{t) contains all fre- 
quency components Ae~*"*. Assuming that the system 
was in a projected random state $ ^^"^ = before the 
perturbation, the wave function after the perturbation 
(t > 0) becomes 



{^E_r\X\^E_r) )) = E iiCCp)) {Ea\X\E,) (42) 

cE^ <Ef 

E (43) 

E^<Ef 

where the statistical relation Eq. ( |30| ) is used. 



* (i) = * ("'(<) + 5* (t), 

* it) = ^ ^°\t) + S4 {t), 



where 



(5* (t) = {-i)e-'^^e(H - Ef)A^Ef , 



and 



<5* (t) = {+i)^Ef Ae{li - Ef)e 



+iHt 



(52) 
(53) 

(54) 
(55) 

(56) 
(57) 
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are the time evolution of unperturbed and perturbed 
vectors. In Eqs. (|55| ) and (p?!), projection operators 
9(H. — Ef) have been introduced to ensure that the ex- 
cited states should be higher than the Fermi energy. 

The linear response of an observable B from all elec- 
trons is calculated as 



dB{t) = 2 Re {s4E,it)B^j^°\t)} 



(58) 



where B = S~^B is the matrix of B in the mixed rep- 
resentation. Then the Fourier transformation of SB{t) 
gives the linear response of the noninteracting many- 
electron system to the perturbation Ae~*'^* , 




5B{t) 



(59) 



where the imaginary part of frequency rj is introduced to 
limit the integration time to a finite value T — — h\5J_n, 
with 5 being the relative numerical accuracy of Eq. (59). 
Here (( ■ )) indicates the statistical average. 



VII. SUMMARY 

We presented a generalized version of the projection 
method for linear and nonlinear response functions devel- 
oped by litaka and others |p^-pT[|. The method can now 
be used with nonorthonormal basis sets such as local ba- 
sis sets, for order- iV total energy calculations. As a result, 
it became possible to calculate the response functions of 
very large systems by applying the projection method to 
the optimized Hamiltonian with a local nonorthonormal 
basis set. 
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